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The spatial structure of the spinor Bose-Einstein conden- 
sates with the spin degrees of freedom is analyzed based on 
the generalized Gross-Pitaevskii equation (GP) in the light 
of the present spin domain experiment on rar = ±1, and 
of the hyperfine state F = 1 of 23 Na atom gases. The 
GP solutions in three- and one-spatial dimensional cases re- 
produce the observed spin domain structures, revealing the 
length scale associated with the existence of the weak inter- 
action of the spin-spin channel, other than the ordinary co- 
herence length related to the density-density channel. The 
obtained domain structure in GP is compared with the re- 
sult in Thomas-Fermi approximation. The former solution is 
found to better describe the observed features than the latter. 

PACS numbers: 03.75.Fi., 05.30. Jp. 



I. INTRODUCTION 

The experiments for Bose-Einstein condensation 
(BEC) in alkali atomic gases, such as 87 Rb 0, 23 Na [§, 
and 7 Li |3| , have been performed under a strong magnetic 
field since its experimental realization in alkali atomic 
gases in 1995. They used the magnetic field to confine a 
BEC system. Because the atom spin direction adiabati- 
cally follows the magnetic field, the spin degrees of free- 
dom are frozen in these magnetic trapping experiments 

Recently, the MIT group has succeeded in creating 
BEC by an optical dipole trap formed by a single in- 
frared laser beam (5]-[7]] in which the spin degrees of free- 
dom are all active (rap = 1,0,-1 of the F = 1 atomic 
hyperfine state for 23 Na). These atoms with the three hy- 
perfine substates simultaneously undergo Bose conden- 
sation, leading to a spinor BEC, a situation analogous 
to superfluid 3 He pj of a neutral Fermion system or a 
triplet superconductor of a charged Fermion system such 
as UPt3 |9|. The spin degrees of freedom play a funda- 
mentally important role for governing their physics. An 
advantage of the present spinor BEC over 3 He or UPt3 is 
that it is a weakly interacting system and we know fairly 
well how the (quasi) particles interact using the knowl- 
edge of atomic physics. These facts allow us to have 
a chance to construct a microscopic many-body theory 



from first principles, using only a few fundamental mat- 
ter parameters. 

Stimulated by the earliest optical trap experiment by 
Stumper-Kurn et al. |^|, general theoretical frameworks 
for describing a spinor BEC were given independently by 
Ohmi and Machida and by Ho [0. They are equiv- 
alent basically. The framework is based on Bogoliubov 
theory which is extended to a vectorial order parameter 
with three components, corresponding to rap — 1,0,-1 
of the F — 1 atomic hyperfine state, giving rise to gen- 
eralized Gross-Pitaevskii (GP) equation. They calculate 
low-lying collective modes such as sound wave, spin wave, 
and their coupled mode and predict various topological 
defect structures, or spin textures. We also note that 
there are many remarkable researches | fl2| with the phase 
separation problem of a two-component BEC consisting 
of different hyperfine states. 

Subsequently, Stenger et al. |jj produce an optically 
trapped spinor BEC with a long cigar shape whose as- 
pect ratio is over 40 and examine whether these three 
component BEC can be either miscible or immiscible. 
The external magnetic field is applied along the long axis 
(the z axis) to see the spin domain formation. After 
releasing the spinor BEC, the Stern-Gerlach separation 
of the cloud is performed to reproduce the original do- 
main structure. Their analyses are based on the Thomas- 
Fermi (TF) approximation, which neglects the kinetic 
energy, of the above generalized GP equations. They 
conclude that the spin-dependent interaction channel of 
the present spinor BEC is antiferromagnetic, rather than 
ferromagnetic in this particular hyperfine state. 

In this paper, we analyze the above experiment by 
Stenger et al. |J in more detail to determine the three- 
dimensional structure of the domain wall including the 
miscible and immiscible spin domain structures. The 
one-dimensional calculation is also done to discuss the 
characteristic lengths. The arrangement of the paper is 
as follows: In next section, we introduce the generalized 
Gross-Pitaevskii equation of a Bose-Einstein condensed 
system with internal spin degrees of freedom. In Sec. 
Ill we simulate the actual experimental system in light 
of the experimental conditions of Ref. |6j and investi- 
gate three-dimensional systems. The properties of a one- 
dimensional system in an idealized situation are explored 
in Sec. [V. The last section is devoted to summary and 
conclusion. 
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1 



II. FORMULATION 



The energy of the system is given by 



The Hamiltonian invariant under spin space rotation 
and gauge transformation is written in terms of the three 
component field operators: ^+1, "to, ^-1, corresponding 
to the sublevels trip = +1,0,-1 of the hyperfine state 
F = 1. Namely, it is given by 

H= fdx^2^l(r)hi(r)^i(T) 

+ yE*!(r)f J t (r)$ 3 (r)f t (r) 

a \ i,j / 



(1) 



where 



fti(r) = - 



n 2 v 2 

2m 



H + V(r)-p(F z ) i;i + q(F*) iti (2) 



is the one-body Hamiltonian, V(r) is the trapping poten- 
tial, i,j — 0, ±1 are spin indices and F a (a = x,y,z) are 
the following 3x3 spin matrices: 



F r 



F„ 



F, 



1 
_.101|, 

v2 V 1 




(3) 
(4) 
(5) 



We have introduced the linear and quadratic Zeeman en- 
ergies §], E ze = E -p(F z ) + q{F%) with p = p + p and 
the Lagrange multiplier po represents the conversation of 
the total spin of the system. 

The interaction constants are related to the two kinds 
of the scattering lengths clq and ai corresponding to the 
total spin zero channel and two channel: 



Airli ao + 2ci2 



9s = 



(6) 
(7) 



Let us introduce the three component order parameters: 
4>i — i^i) with i = x, y, z. Following the standard pro- 
cedure p3| , we can derive the Gross-Pitaevskii equation 
for describing 4>i(r) from the Hamiltonian Eq. ([j]): 



(r) 



a j.k,l 



(8) 



E - fdrJ2{ri(r)hi(T)<Pi(r)} + jdv{E n {v) + E s (r)}. 



where 



E n (r) = fJ2\Mr)\\ 

i 

a \ i.j 



(9) 

(10) 
(11) 



We will use the following notation in Sec. |Vf (i(r) = 



where n(r) = ^ 



III. THREE DIMENSIONAL CASE 

In this section we analyze the three-dimensional sys- 
tems in comparison with the actual experiments done 
by Stenger et al. ||. The main interests are (1) the z- 
position of the domain wall, and (2) the radial shape 
including the mutual overlapping of the spin domains. 



A. Experiment 

In order to establish a model system we briefly explain 
their experimental conditions. The condensate is under 
the magnetic field B(z) which is applied along the z axis. 
Since the atoms are in the hyperfine F = 1 state, the con- 
densate consists of the three components: <p + i,(pQ, and 
<f)—i. The non-uniform magnetic field B(z) is character- 
ized by its field gradient: p' oc dB(z) / dz. The p in E zc 
is set to zero at z = and p = p'z. The stronger-field 
side of the cigar-shaped condensate is filled with the +1 
condensate. As the system follows the total spin conver- 
sation, the opposite side is filled with the —1 condensate. 
Stenger et al. Q observe the column density distribution 
of the three spin components and estimate the interac- 
tion constant g s from the z-position of the domain wall 
between the spin components. The result depends not 
only on p, but also on the coefficient q. q cx Bq where 
Bq ~ B(z) is the base magnetic field. The experiments 
are done for various p's and g's. 

They estimate the relationship of p and q with the z 
position of the domain wall by approximating the cigar- 
shaped system as a one-dimensional one with uniform 
density, which is assumed to be 2/3 times the peak den- 
sity. The kinetic term is also ignored. This allows us 
to draw lines in p-q phase diagram as a function of the 
interaction constant g s (in their paper, C2). 
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B. Calculation 

In order to check whether or not their analysis based 
on the above assumptions (1-dimensional TF approxi- 
mation, ldTF) agrees with the full 3-dimensional calcu- 
lations of GP equation, we determine the spatial pro- 
files of the rotationally-symmetric (around the z axis) 
spinor condensate by solving the GP equation in the 
r-z space using the difference between the scattering 
lengths ai — ao = 0.19nm ( oc g s ) and the peak density 
"peak = 4.35 x 10 20 m~ 3 they estimated. Assuming this, 
we compare the following three calculations by ldTF, 
by 3-dimensional TF approximation (3dTF), and by 3- 
dimensional GP equation (3dGP). In ldTF, it is assumed 
that the density is |fT. p0 ak- 

We first discuss the domain structure along the r di- 
rection and the mutual overlapping between the three 
components. One example (p' = 1.0Hz//ito, q = 2.0Hz) 
of the calculated results with 3dGP is depicted in Figs. 
H| and ||(a). There is a large overlapping region consist- 
ing of the +1 and the -1 components, between the re- 
gion consisting of the +1 (or -1) component only and of 
the component only. In other words, the double peak 
structures of the —1 component and +1 component are 
situated both at the z > side and z < side in Figs. 
|l| and ||(a). This feature accords with the observation 
(Fig. 3 in ||). Fig. ^(a) also shows that the condensate 
itself allows the overlapping of 3 components at z = 
in certain values of p' and q. This existence of three 
components at z = does not occur in calculations with 
3dTF [e.g. Fig. g(b)]. Usually the difference between the 
results of TF approximation and of full GP equation is 
considered to be an order of the coherence length ~ 1/jm 
here, but the difference between the two treatments is as 
long as 50/im, which is far longer than expected. This 
behavior also occurs in the one-dimensional system and 
is discussed thoroughly at Sec. IV. (This overlapping at 
r = does not occur at larger q, but the length scale 
stays as long as 10^m.) 

Another interesting difference between the results from 
3dTF and 3dGP exists at the crossing point between the 
component and the +1 component. As seen in Fig. 
||(a) , the z coordinates of the crossing points do not move 
appreciably between the two cases: r = and r = 2. 
However, the crossing point moves significantly in TF as 
seen from Fig. |^(b). The trace of these crossing points in 
the r-z plane for various radial harmonic potentials are 
depicted in Fig. ||. 



C. Phase diagram 

To compare the z-position of the domain wall between 
the component and the 1 component, we integrate the 
density profile e.g. Fig. |l| along the r surface. Figure |] 
(a) shows an example of the results in 3dGP and Fig. 
[| (b) shows that in 3dTF. When the density profile is 



integrated along the r direction this produces further 
overlapping between the components especially in 3dTF 
jcompare Fig. |^(b) with Fig. ^(b)]. As seen from Fig. 

(a) , the z coordinate of the crossing point in 3dGP is 
slightly smaller than the 3dTF case in Fig. ^ (b). 

The acquired z value of 3dTF for various p 1 and q val- 
ues are shown in Fig. §(a). The ldTF line p = 2^/MT^ 
is also plotted. The value "50.7" is derived uniquely 
from the assumed scattering lengths and the peak den- 
sity (x2/3). No significant difference is seen. The lower 
branch lines in Fig. ^ comes from the density profile along 
the z axis; When z is large, the value "50.7" (oc density) 
must become much smaller. Figure |s](b) compares the 
ldTF and the 3dTF. 

The differences in p-q diagram between ldTF, 3dTF, 
and 3dGP are small. Therefore we conclude that the 
estimate in Ref. is correct. But the shapes of domain 
walls are not simple as shown in Figs. and ||. 



IV. ONE DIMENSIONAL CASE 

We have seen in the previous section that the difference 
between the density distributions by the TF result and 
GP result is much wider than the usual coherence length 
£ < lpm. In this section, we take up the one-dimensional 
system to investigate why this is so, taking into account 
the effects of the kinetic term and the interaction g s term. 
The magnetic field coefficients (p and q) and the optical 
potential V(r) are set to zero in this section. The kinetic 
term, the interaction (g n and g s ) terms, and the chemical 
potential are retained. 



A. Positive g s 

Figures ^(a) and (b) show the cases when g s is posi- 
tive. The ratio of the components is fixed at each end of 
the system; (C+i, Co, C-i) = (0, 1, 0) at the left hand edge 
and (1,0,0) at the right hand edge of the system. We 
impose an additional condition that the C-i component 
is zero in Fig. ^ (a). In this case, the overlapping between 
Co and C+i is seen to be an order of a few micrometers. 

This length scale can be explained by £ s e 



Sirn(a2 — ao) ' 

which is about \.2p,m in the present case. Figure | (b) 
shows that when no component is suppressed, the over- 
lapping region becomes wider and is subdivided into the 
wider region (— 10/ito < z < 5/im) where £-i = C+i and 
the narrow region (5fj,m < z < 10/zm) where C-i 7^ C+i- 
The latter region is characterized by the coherence length 

The difference between Figs. |y(a) and (b) is explained 
by considering the E s term of the total energy [Eq. (g)]. 
The relative phases of £'s (equivalently, of </>'s) are deter- 
mined under the condition that E s is minimized as shown 
in Appendix |X[ We can write E s as 
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a \ i,j J 

= ^[-(IC+il ± IC-il) 4 + 2(|C+i| ± IC-il) 2 ]. (12) 

The upper and lower signs correspond to the positive and 
negative g s respectively. Equation ( |l2| ) shows that E s is 
minimum on the line |C+i| = |C-i| when g s > 0, and on 
IC+i| + |C_i| = lwhen. 9s <0. 

Figure ||(c) shows the landscape of E s when g s > 0. 
The gray lines (a) and (b) in Fig. ^|(c) correspond to 
Figs. g(a) and (b) respectively. The line (b) goes along 
the minimum energy line |C+i| — |C-i| &t nrs t an d this 
corresponds to the z = — 10/im to 5/im region in Fig. 
^(b). Both the lines (a) and (b) climb up the hill in the 
E s landscape from E s = to E s = n 2 g s /2. This process 
is controlled by the kinetic term and the interaction g s 
term in the GP equation Eq. (g), and the length scale is 
an order of £ s . 

The spatial variation in Fig. |^(b) is also understood by 
so-called d- vector [||, which expresses the spin structure 
of superfluid. At z = — 10/im and z = 10/im the d- 
vectors are given by d — z and d — x + iy respectively. In 
the region from z — —10/im to 5jum the d- vector rotates 
from the z-direction to the x-direction. Since these states 
are energetically degenerate, this length scale is governed 
by the system size. On the other hand, in the region 
from z — 5/im to z = 10/im the d-vector is described 
by d = x + ity where t varies from t = {z = 5/im) to 
t = 1 (z = 10/im). This state change produces the energy 
variation related to the interaction g s , thus is governed 
by the spin coherence length £ s . 

B. Negative g s 

Figure H shows the cases when g s is negative. In 
Fig. (a), the ratio of the components is fixed to 
(|C+i],|Co|,|C-i|) = (0,0,1) at the left hand edge and 
(1,0,0) at the right hand edge of the system. As shown 
by the gray line (a) in Fig. |?| (c), the system goes along the 
minimum E s line |C+i| + |C-i| = 1- Therefore the length 
scale is not controlled by £ s and it becomes as long as 
the boundary conditions allow, that is, the length scale 
is determined by the boundary condition. This gentle 
shape contrasts remarkably with Fig. ^ (b) which shows 
the system when the ratio of the components is fixed to 

(|C+i|, |Co|, | C x 1) = (0,0, 1) at the left hand edge of the 

system. The quick change around z = —10/im to — 7/xm 
is controlled by £ s in a similar reason of the previous 
g s > case. It is interesting to note that the £_x compo- 
nent spontaneously appears to minimize E s \ the system 
without restriction does not go along the |C-i| = line. 



C. Length scales 

The length scale of the density variation in the sys- 
tem without the spin freedom is determined by £ = 
y/l/{%ima) where a(> 0) is the s-wave scattering length. 

This is equivalent to £„ = Svn ^l 2+ ^ j(. 2a 2 + a o > 0) in 
the system treated in this paper. The origin of £ n is the 
competition between the kinetic term and the interaction 
term with g n in the GP equation. 

There are three types of the length scale in this system; 

(1) £ n . This becomes dominant when the total density 
varies. This is not significant in these one-dimensional 
systems because their densities are almost uniform. In 
three-dimensional calculations, the total density profile 
along the r axis in Fig. [l], which is more gentle than 
that in 3dTF (no corresponding figure in this paper), is 
explained by £„. 

(2) £ s . When the total density does not vary signifi- 
cantly, this length scale becomes important. The origin 
is the competition between the g s term and the kinetic 
term in the GP equation. For example, Fig. ^](a) and the 
left hand side of Fig. @(b). 

(3) The length controlled by neither £ n , nor £ s . When 
the spin state of the system varies along the lowest energy 
line of E B , the length scale of the density variation is 
governed neither by £„, nor by £ s . The characteristic 
length is controlled by the boundary conditions in the 
one-dimensional cases: The left hand side of Fig. |(b) 
and Fig. 0(a). They correspond to the line (b) in Fig. 
^(c) and the line (a) in Fig. 0(c) . 

As for the previous 3-dimensional calculations, the 
type (3) of these length scales explains the difference be- 
tween the results of the 3dGP and the 3dTF in Figs. |(a) 
and (b). The magnetic field parameters p and q deter- 
mine the characteristic length when 3dGP is used. It is 
neither £ n nor £ s . We note that the (+1 and £_i compo- 
nents have the finite density at z = in Fig. ||(a). This 
long length scale is explained by the minimum energy 
line of the g s term in the GP equation (or E s ). This is 
the |C-i| = |C+i| line of Fig. §(c). 

V. SUMMARY AND CONCLUSION 

We have investigated the spinor BEC systems based on 
the generalized GP equation extended to the cases where 
the BEC has the spin degrees of freedom. The spatial 
structure of the domain wall which is the interface of the 
different spin states of BEC is analyzed by solving the GP 
equation both for rotationally-symmetric (around the z 
axis) three-dimensional and one-dimensional cases. The 
former case simulates the actual experimental situation 
of Stenger et al. || . Our calculations show that the sim- 
ple TF approximation taken by Stenger et al. can be 
justified and yields the correct value of the interaction 
constant of the spin channel which turns out to be an- 
tiferromagnctic in the present hyperfine state F = 1 of 
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23 Na atoms. The full GP solutions both in three- and 
one-dimensional cases reveal a long length scale associ- 
ated with the interaction of the spin channel, yielding the 
large overlapping region between the immiscible compo- 
nents (e.g. and +1). 
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APPENDIX A: RELATIVE PHASES OF <j> 

The relative phases of the three component condensate 
wavefunctions are determined so that the energy density 
of the g s term 



a \ i.j / 



(Al) 



is minimized under the condition that the amplitude of 
(f>'s are fixed. This is because the other terms in E are 
not affected by the choice of the phases. We assume that 



,</»o,</>-i) = (/3ie J71 ,/5 ,/?-ie^- 1 ) 



(A2) 



where /3j(> 0) and 7, are real numbers. The amplitude 
Pi ar e fi xed and we determine 7, to minimize E s . From 
Eq. (Al) we obtain 
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+ (l^i| 2 -|^l)| 2 

= ^[(201 cos 7l + 2/?_ 1 cos 7 _ 1 ) 2 

+ (2/3! sin 7i-2/3_i sin 7 _i) 2 ] 
+(/3 2 -/3 2 x) 2 
= 2$\ffi + !3\ + 2/3i/3_i cos( 7l + 7 _ x )] 
+(/3 2 -/3 2 1 ) 2 . (A3) 

When g s > 0, E s is minimized with 71+7-1 = n. There- 
fore, we can take 6 as 



(<j)i,<f> ,<f>-i) = (/3i,/3 , ~/3-i). 
When g s < 0, 7 i + 7 _i = 0. Therefore, 



(A4) 
(A5) 
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FIG. 1. The density plot of the r-z plane. The dashed 
lines, the solid lines around 2 = 0, and the dotted line mean 
the densities of the +1, 0, and —1 components. The upper 
solid curve at r = indicates the total density at r = 0. 
The peak density is 4.35 x 10 20 m~ 3 , p' = 1.0Hz/^m, and 
q — 2.0Hz. The optical confinement is assumed as the har- 
monic potential given by V(r) — ^{(2Tvv r r) 2 + {2-kv z z) 2 }, 
where v z = 15Hz and v T — 900Hz. The total particle number 
is 8.71 x 10 5 . 
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FIG. 2. The density profiles of the condensate. The peak 
density, p, q, and the shape of the harmonic confinement 
are same as in Fig. |l| except for the total particle number 
8.99 x 10 5 in (b). The dashed lines, the solid lines around 
z = 0, and the dotted line mean the densities of the +1, 0, 
and —1 components. The upper solid curves indicate the total 
density, (a) The density profile at r = and 2fim. The ±1 
components exist even at z = 0. We call the crossing between 
the and +1 components (circles) the domain wall. The z 
coordinates of the domain wall are z = 19.5/im at r = and 
z — 19.2/im at r = 2/tm. (b) The density profile when cal- 
culated with 3dTF. As r increases, the domain wall (circles) 
between the component and ±1 components moves toward 
z — 0. The z coordinates of the domain wall are z — 24.1pm 
at r — and z = 16.2/xm at r = 2pm. 




0.08 0.12 0.16 

z/Ztf 

FIG. 3. The cross-lines between the component and the 
+ 1 component for v r = 30, 75, and 900Hz. The solid line, 
which is independent of v r , is the line when TF approxima- 
tion is used. The 900Hz line correspond to Fig. ^. Each axis is 
normalized by the TF radius and length, which are defined as 
(peak density) x g n = f(2-KV r R TF ) 2 = f (2tt^Z T f) 2 . The 
Rtf is 82.2, 32.9, and 2.74/tm for v T = 30, 75, and 900Hz re- 
spectively. Ztf is 164/im. 
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(a) The r-integrated density profiles for p — 1.0Hz//im and q = 2.0Hz. The overlapping of these integrated 
densities simply reflect the shape of Fig. [| (a). The z coordinate of the domain wall is 19.3/^m. (b) The corresponding result 
for 3dTF. The overlapping of these integrated density comes from the r-dependent shift of the domain wall depicted in Fig. ^. 
The z coordinate of the domain wall is 2l.5(im. These figures should be compared with Fig. 3 in |H|. 
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FIG. 5. The p-q diagram, (a) The comparison between the 
the results in ldTF and 3dGP. (b) The comparison between 
the results in ldTF and 3dTF. The magnetic field parameters 
are p' = 0.2, 0.5, 1, 2, and 4Hz//im, and q = 1, 2, 5, 10, 20, and 
40Hz (not all of the combinations are used). The ldTF line 
is p = 2y/507Tq. 
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FIG. 6. The one-dimensional systems with positive 
g 3 (ao = 2.75nm, = 2.94nm). The chemical po- 
tential /i is 4.35 x 10 20 m" 3 x g n . (a) The density vari- 
ation of the three components when the boundary con- 
dition (|C+i|,|Co|, |C-i|) = (1,0,0) at z = -10/xm and 
(0,1,0) at z = 10/im is imposed. We suppress f_i over 
whole range of z. (b) The density variation of the three 
components when the boundary condition is imposed that 
(iC+ij, ] Co 1 , |<-i|) = (1,0,0) at z = -lOfim and (0,1,0) at 
z = 10/^m. (c) The landscape of E B . The labels (a) and (b) 
correspond to the above figures (a) and (b) respectively. 
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FIG. 7. The one-dimensional systems with negative g s 
(do = 2.75nm, a2 = 2.56nm). The chemical potential fi 
is 4.35 x 10 20 m -3 x g n . (a) The density variation of the 
three components when the boundary condition is imposed 
that (|C+i|, |Co|, |C-i|) = (0,0,1) at z = -10/xm and (1,0,0) 
at z = Wfim. (b) The density variation of the three 
components when the boundary condition is imposed that 
(|C+i|, ICo], IC-il) = (0,1,0) at z = -Wfim and (1,0,0) at 
z = 10/im. (c) The landscape of E a . The labels (a) and (b) 
correspond to the above figures (a) and (b) respectively. 
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